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Abstract 

We investigate a simplified model of two fully connected magnetic systems maintained at different 
temperatures by virtue of being connected to two independent thermal baths while simultaneously 
being inter-connected with each other. Using generating functional analysis, commonly used in 
statistical mechanics, we find exactly soluble expressions for their individual magnetisations that 
define a two-dimensional non-linear map, the equations of which have the same form as those 
obtained for densely connected equilibrium systems. Steady states correspond to the fixed points 
of this map, separating the parameter space into a rich set of non-equilibrium phases that we analyse 
in asymptotically high and low (non-equilibrium) temperature limits. The theoretical formalism is 
shown to subvert to the classical non-equilibrium steady state problem for two interacting systems 
with a non-zero heat transfer between them that catalyses a phase transition between ambient 
non-equilibrium states. 
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I. INTRODUCTION 



The study of energy transport between interacting non-equilibrium systems, each in con 



tact with heat baths maintained at different temperatures, has a long history [l-|5[ span- 
ning a multitude of non-equilibrium systems. Two macroscopic systems in thermal contact 
with each other will exchange energy until equilibrium is reached as long as they are kept 
thermally isolated from the environment. However, if systems are coupled to independent 
thermal reservoirs that are maintained at two different temperatures, they can be kept 
indefinitely in a non-equilibrium steady state (NESS) without any perceptible change in 
macroscopic local properties. While the definition of a 'non-equilibrium temperature' has 
remained an open question for long, not the least due to the absence of an inherent Hamil- 
tonian description, generalised manifestations of the fluctuation-dissipation theorem (s 9 1 
have often addressed the issue of non-zero thermal flux across non-equilibrium quasistatic 
systems that are driven by external time-dependent forces. The common topic of interest in 
all these studies has been the depiction of the temperature profile as a function of the free 
energy differences of connected systems and the dependence of associated heat flows on this 
'non-equilibrium temperature'. 

A recent series of studies have specifically focused on stochastic fluctuations (2- 4 ] as 
the initiators of external forcing that are simultaneously capable of maintaining the non- 
equilibrium energy traffic across the connected (non-equilibrium) systems, both in lin- 
ear 2|, |5] and non-linear regimes 10|. The models considered include harmonic chains (in 
one, two and three dimensions) and harmonic crystals with (mostly) white-noise Langevin 
reservoirs connected at the ends of the chains/ crystals. All these works indicate an analyt- 
ical dependence of the non-equilibrium heat transfer on the dimensionality of the system, 
categorically characterising the phase space of the combined system into either of the three 
standard universality classes - ballistic, diffusive and anomalous (sub or super- diffusive). The 
non-equilibrium manifolds inclusive in these studies are also known to abide the celebrated 



'additivity principle' 



111. 



Interestingly enough, none of these studies considered externally forced spin systems as 
the perturbed basis, be it in the form of the paradigmatic spin-1/2 Heisenberg system, or 
its equivalent spherical n-spin model. In view of the wide array of spin and spin-glass type 
systems undergoing non-equilibrium energy exchanges that lead to steady state that exhibits 



metal-insulator-like conduction, this appears a unique exclusion. One must mention [12| as 
an exception in this regard; where starting from spinless quantum dots, the authors were able 
to prove the existence of two-particle scattering states leading to a steady non-zero Landauer 
current. While of a different flavour, another work worth mentioning is the interpretation 



of the replica method as a two temperature system that exhibits different dynamics re 



ated 



to the corresponding effective temperatures of annealed and quenched variables |13Hl5|. 

Theoretical implications aside, the two-temperature formalism we examine here has a 
multitude of direct applications in condensed matter systems. Many solids, including the 
high-T c cuprates and pnictides superconductors, have layered structures where intra- and 
inter-layer couplings are different, leading to the existence of two different temperature baths 



in an inter-connected spin system [18|. As would be detailed later, a corollary of our model 
addresses the simplified case of a temperature gradient perpendicular to the layers that 
effectively embeds anharmonicity in the interacting spin chains. Identically, superlattices 
subjected to anharmonic (and often aperiodic) forcing while being maintained at two differ- 



ent temperatures, also fall under the more general purview of the model presented here 



19|. 



Interestingly, other potential applications in the study of dark matter, where it has been 



suggested 



161 ] that normal and dark matter are coupled to each other through gravity while 



remaining in contact to thermal reservoirs at different temperatures in their respective sec- 



tors. Some 12| even propose that the actual equilibrium temperature of each sector can be 
different. 

As a major departure from the existing trend, in this article our NESS connected sys- 
tems are deterministically forced and time independent to start with. More specifically, 
we consider two mutually connected spin-1/2 chains that are individually attached to two 
separate heat reservoirs that are maintained at two separate fixed temperatures through 
steady supplies of external energy. Our interest here is in the magnetisation regime of 
each system under the action of mutual (heat) flux exchange and how such non-equilibrium 
magnetisation connects to equilibrium steady states, a direct allusion to the alternative 
fluctuation-dissipation formalism. 

The experimental analogue of our theoretical framework is that of two inter-connected 
magnetic systems adiabatically maintained at two different constant temperatures by virtue 
of being connected to independent temperature baths. The two magnetic systems that are 
suffused with non-zero magnetisations at the initial time instant t = are then allowed to 



interact via magnetic couplings and left to evolve in time. Physically, this can be accom- 
plished by bringing two such magnetic systems close together in vacuum or by connecting 
them through some thermal insulating material. The steady states of this composite system 
will then be analysed under synchronous discrete Markovian dynamics. 

The microscopic model and its dynamical rule will be presented in section [III the solution 
of which will be a coupled two-dimensional dynamical map for the magnetisations of the 
two subsystems. As we are interested mainly in the steady state solutions, we will analyse 
in some depth the fixed points of this map on section IIII1 Section IIVI contains a stability 
analysis of the global paramagnetic phase of the model. Further discussion of the results, 
their relevance and future directions are presented in section [V] 

II. THE DYNAMICAL MODEL 

We start by considering a two-dimensional fully connected Heisenberg lattice where each 
node represents a spin-1/2 particle. We then identify two different sets of nodes, a and r, 
which we refer to as the a and r components of the system (subsystems). We ascribe the rule 
that nodes in each one of the components are fully connected to their own separate thermal 
reservoirs and no particle is connected to more than one particle of the other subsystem. 
The reservoirs to which the a and r components are coupled to temperatures T a and T T , 
respectively. We assume that the cx-component has N a nodes <7j G {±1} and, similarly, the 
r-component has N T nodes Tj G {±1}- 

The system dynamics will be taken as discrete and defined by the local stochastic update 
rules for each individual spin-node by the probabilities 



with P a = 1/T a , a = <j,t, where we are using units such that the Boltzmann constant is 
ks = 1. The local fields acting on each spin are defined as 



P[<n{t + 1)] oc exp {/3 a ai(t + l)hf[a(t), r(t)]}, 



(1) 
(2) 



P[Ti(t + 1)] cx exp {p rTi (t + l)hj[a(t), r(t)]} 




T ._/ . 



(3) 



(4) 
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FIG. 1. Visualisation of the model as two planes a and r. The long-ranged intraplane interaction 
is J and the short-ranged interplane interaction is J' . 



The above description suggests that constituents of each component interact with one an- 
other via long-range (mean-field) interactions J a and with spins from the other compo- 
nent locally via the short-ranged (local) interactions J' ab . For simplicity, we assume that 
Jcr — Jt — J, J'ot = J'to = J' an d that the number of spins in both components is the same 
N a = N T = N. 

A helpful way of visualising this system is depicted in Fig. [T]where we consider the a and r 
components lying on two different parallel planes. We will use this graphic interpretation as 
it conveniently separates the global and local interactions and simplifies the nomenclature we 
will be using. The long-range interaction defined by the coupling J will be termed intraplane 
interaction, while J' will define the local and symmetric interplane interaction. 

The system is updated synchronously, meaning that at each time step all the spins are 
updated simultaneously; the generating functional method 20| then allows for the dynamics 
to be solved exactly. This results in a two-dimensional non-linear dynamical map for the 
magnetisations of the different components of the system given by 



nf{t + l) = (tanh/3 (T /i ,7 (t)) T(i) , 
m T (t + 1) = (tanh/3 T /i T (t)) CT(t) , 



(5) 
(6) 
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where 



h a (t) = Jm a (t) + J'r(t), (7) 
h T (t) = Jm T (t) + J'a(t), (8) 

and where averages (•) are with respect to the corresponding probability distributions 

P[a{t)] = \[l + a{t)nf{t)], (9) 
P[r(t)]= 1 -[l + r(t)m^t)}. (10) 

The steady state solutions are given by the fixed points of this map 

m a = (tanh/3 ff (Jm a + J'r)) T , (11) 
m T = (tanh (3 T (Jm T + JV)) CT , (12) 

which are exactly the same equations that describe, in equilibrium, the magnetisation of a 
spins system with long-range interactions subjected to a random magnetic field. These fixed 
points are analysed in the following section. 



III. STEADY STATES 



Steady state solutions of the dynamics of our system will be given by the fixed points of 
the two-dimensional map defined by equations (JSJ) and ([6]). A notational simplification of 
those equations can be achieved by using the reduced variables (3 a = J'(3 a and J = J/J', 
resulting in the following map 

r = ^^(l + sm r )tanh/3 ff (ym a + /), (13) 

s 

f T = - (1 + sm a ) tanh f3 T (jm T + s) , 

s 

with s G {±1}- The dynamics now depends only on three parameters j3 a , j3 T and J and the 
initial conditions on the magnetisations. Note, though, that (3 a and /3 T may assume negative 
values depending on the sign of the interaction J' and cannot be interpreted straightfor- 
wardly as simply rescaled inverse temperatures. 
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It is possible to obtain the fixed points of these equations numerically for any parameter 
set. Still, it is interesting to study the behaviour of the system in the limits of high and low 
temperatures, which can be obtained analytically and exposes a very rich phase diagram. 
The latter represents different types of interactions between the systems components and 
provides insights into the corresponding macroscopic properties. 

The system, as defined, is symmetric under exchange of component labels cx -H- r. There- 
fore, without loss of generality, we choose the r-component to be connected to a thermal 
reservoir at high temperature /3 T — > 0. Then 



resulting in the fixed point 



f « ~$ T {jm T + m CT ), (14) 



m T = (3 T m a , (15) 



at leading order in /3 T . 

Naturally, at high temperature the intraplane interaction J is not sufficiently strong to 
influence the spins alignment and disappears from the equation. The interplane interaction 
still survives, although weakened by the thermal disorder, linking the magnetisations of 
the two components via a proportionality constant f3 T ; the relative sign between the two is 
determined by the nature of the interaction J' (ferromagnetic or antiferromagnetic) in an 
obvious way. 

When the r-component temperature becomes infinite, with /3 T being exactly zero, the 
intraplane interaction becomes too weak to allow for the interaction with the cx-component 
to affect the r-component alignment, which then leads to a paramagnetic phase with mag- 
netisation m T = 0. It is easy to see that consequently the cx-component becomes equivalent 



to a mean-field Ising model in a random field of zero mean [22 1. 

Keeping the r-component's temperature fixed at a high but finite value, one can re- 
examine the system's properties at various temperature limits of the cx-component. For 
the case of high temperature (low one obtains a proportionality relation between the 
magnetisations 

m a = a m T m T = f3 T (3 a m T . (16) 

As we are assuming both temperatures to be high (both /3's small), the only solution is 
that both magnetisations vanish at high enough finite temperatures. This signals a transition 
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FIG. 2. Fixed points of the magnetisation m a for J = 2 and $ a = 3 

to what we conventionally refer to as the global paramagnetic (GP) phase. The GP phase 
is a ubiquitous fixed point, present for any temperature, but its uniqueness and stability 
are temperature-dependent. In this section we will focus on its uniqueness property, i.e., 
analyse the existence of other possible phases in these temperature limits. The issue of its 
stability will be addressed mainly in the next section. 

One notable feature of the high-temperature regime in r is that it is rich enough to develop 
a total of either three or five fixed points as the a-component's temperature decreases, 
depending on the values of the intraplane interaction J'. This can be obtained from the 
relation between f a and m a of Eq. ( {TBI by utilising the high temperature relation ( fT5l) and 
the intersection of f a with the line m a . Figure [2] shows an exemplar plot for J = 2 and 
$ a — 3. The graph shows five fixed points: a stable solution at m" = 0, two stable solutions 
at high magnetisation values (largest) and two unstable solutions at intermediate values. 

For r at high temperature and at the limit \(3 a \ — > oo one obtains 

rrf = (sgnf3 a (^Jm a + rj^ = sgn J' • <^sgn (jrrf + r) ^ . (17) 

In this case, one has to consider the various possible values for the intraplane interaction 
J. For \J\ < 1 

rrf = m T sgn f3 a = m T sgn J', (18) 

which, combined with the proportionality between m T and m cr , gives rise to a single solution 
m T = m" = indicating that the intraplane interaction is too small to align the spins of the 
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er-system. A straightforward, but lengthy consideration of the possibilities for | J\ — 1 gives 
the same result. This means that there are situations where the GP phase is the only phase 
present at the high-temperature /3 T — > limit for both low and high T CT ; which reveals the 
importance of the relative (non-dimenionalised) interaction strength J in defining the state 
of the combined system, rather than their actual individual values. 
Non-GP solutions do exist for \J\ > 1. When |Jm CT | > 1 we have 

m" = sgn J ■ sgn m" , (19) 

and m a becomes independent of m T . Ignoring the GP solution, we are left with m cr = 
±sgn m" depending on the sign of the intraplane interaction. 

It is convenient to consider several different cases: (i) A positive sign gives m a = ±1, 
meaning that the intraplane interaction is sufficient to align the a component spins com- 
pletely, independently of the heat flow from the r to the cr-components. (ii) A negative 
sign results in oscillations between ±1, a cycle of period 2 for m cr instead of a fixed point 
as depicted in Fig [3j This occurs, for instance, in the case antiferromagnetic intraplane 
interactions and is characteristic of synchronous update dynamics, (iii) Finally, the case 
Jm cr = ±1 has as the non-GP solutions 

m° = ± ^i, (20) 
2-Pt 

which exist when the interactions obey the temperature-dependent condition J = (2 — 
T ) sgnJ'. This means that for each /3 r value only one solution exists on the J line. 

As the equations are completely symmetric with respect to exchanges of a and r, the 
only remaining case is the one where both systems are at low temperatures with T T ,T a 
1, or equivalently |/3 r |, \/3 a \ 3> 1. Once more we consider separately cases with different 
magnitudes of the interactions: 

(a) For | J\ < 1 we have 

/ CT =m r sgnJ', (21) 
f T = m a sgnJ'. (22) 

In this limit, all points of the (m CT , m r )-plane are either period-2 cycles or fixed points; the 
fixed points correspond to the line m a = m r sgn J'. 
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FIG. 3. Dynamics of the 2D map of the magnetisations for f} T = 10 2 , /3 CT = 10 2 , J = —5 and 
J' = 1CT 2 . The period 2 cycle in rricr sets in already at the first iteration, while m T remains at zero 
mainly due to its very high temperature. 

(b) When \J\ = 1, the same solution is obtained for \m a \, \m T \ < 1, but when \m a \ = 
\m T \ = 1 the fixed points m CT = m T sgn J' occur only for sgn J' = J; otherwise the solutions 
are period-2 cycles. 

(c) The | J\ > 1 case is more complicated to analyse. Here the intraplane interaction is 
relatively strong and many different solutions can be found. We have to take into consider- 
ation one of the three cases for each of the magnetisations (a, b — a, r): 

| Jm a | > 1 =^> f a = sgn J ■ sgn m a , 



\Jm a \ 



r 



(23) 



2sgn J' - J 
\J m a \ < 1 =>. f a = m 6 sgn J'. 

(i) The case of |Jm°"| > 1 and |Jm r | > 1 results in two independent equations of the 
form (fT9|) for each one of the magnetisations, giving completely aligned states or period-2 
cycles for all points with \m a \ > 1/\J\ and \m T \ > 1/\J\- (ii) Another simple case is when 
iJm""! < 1 and |Jm r | < 1, or equivalently \m a \ < 1/\J\ and \m T \ < 1/\J\, which reduces to 
the same set of solutions as for | J\ < 1. 

(hi) When |Jm CT | = |Jm r | = 1, or \m a \ = \m T \ = 1/\J\, the equations become 



m 



2sgn J' - J 



m 



2sgn J' - J 



(24) 
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FIG. 4. Partial phase diagram for the limiting cases of high and low temperatures. The figure shows 
only the first quadrant as the other three are simply its symmetrical images reflected about the 
axes. On the dashed axis, the corresponding magnetisation is zero as the temperature is infinite. 

which exists when the interactions obey the relation 2sgn J' — J = ±1. For the plus sign, 
we have fixed points where m" = m T , while the minus sign gives m" = —m T . 

It remains to analyse the mixed cases with: (iv) |Jm a | > 1, \ Jm b \ < 1; (v) |Jm a | > 1, 
\Jm b \ = 1; (vi) \Jm a \ < 1, \Jm b \ = 1. 

Cases (iv) and (v) can only have non-GP phases for m a = ±1 and sgn J = +1 and a 
2-period oscillatory behaviour when sgn J = —1 with m a oscillating between +1 and — 1. 
In both cases, m b will be a multiple of m a according to the corresponding value given by 
one of the equations ( l23i) . Finally, all points in case (vi) are either fixed points, when 
2 — Jsgn J' = +1, or 2-period cycles, when 2 — Jsgn J' = — 1 

Figure HI together with table [J presents a summary of the possible asymptotic phases 
of the system in the limits we analysed. Intermediate regimes are more difficult to analyse, 
but the equations are simple enough to be solved numerically for any parameter values. 

Due to the non-linearity of the equations, we carried out extensive simulations with many 
different sets of parameters looking for an indication of chaotic behaviour, but did not find 
any. Presumably this is due to the fact that the equations of macroscopic order parameters 
result from a stochastic process fluctuating very closely to the typical values at the large 
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TABLE I. Phases of FigH 



Phase 



Description 



T a ,T T > 1 



Global Paramagnetic (GP) only (m CT = m T = 0) 



m 



II 

T a > 1 
a,b = a,T 



\J 



f3 a m b 
< 1 



\J\ > 1 



GP-only, 
\Jm b \ < 1 
Jm b = ±1 



< 



\Jm b \ > 1 



GP-only 

m b = ±sgnJ'/(2-^ a ) 
J > =>• |m 6 | = 1 
J < =>• Period-2 Cycle in m b , 



\J\ < 1 



\J\ = 1 



\m 



\m 



1 



III 

T a ,T T <C 1 



|J|>1 



m a = ?n T sgn J' =>■ Fixed Points 
all other points =>■ Period-2 Cycles, 
\m a \, \m T \ < 1 =4> m' = m r sgn J' 

J = sgn J' =>■ m' = m T sgn J', 
J = —sgn J' =4- Period-2 Cycles, 
J > =>- Totally Aligned, 
J < Period-2 Cycles, 
Same as | J| < 1, 

A = 2 sgn J' — J = ±1 =>m" = Am T , 

All other cases =>• GP-only, 
J > ^ m a = ±1, ?n b = ?n a sgn J', 
J < => Period-2 Cycles, 
J>0 =4> m a = ±1, m b = m a /(2sgn J' - J), 
J < Period-2 Cycles, 
2 — J sgn J' = +1 All fixed points, 
2 — J sgn J' = — 1 =>• All period-2 cycles. 



| Jm '!, | Jm "! > 1 

| Jm a \, \ Jm a \ < 1 
|J?n CT |, |Jm CT | = 1 

| Jm a \ > 1, | J?n b | < 1 

| Jm a \ > 1, | Jm b \ = 1 

I Jm a \ < 1, | Jm b | = 1 



system limit; as such the system is not amenable to chaotic behaviour. 
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IV. STABILITY OF THE GP PHASE 



As previously noted, there are two different types of paramagnetic phases in the presented 
model: (i) the GP phase where m" = m T = 0, and the (ii) a-component paramagnetic phases, 
where m a = 0, with a either a or r and the magnetisation of the other component being 
non-zero. 

The GP phase is ubiquitous in the studied model, being always a solution for any combi- 
nation of parameters. However, depending on the values of these parameters, the stability 
of this phase can drastically vary; it can be analysed by calculating the eigenvalues of the 
Jacobian matrix of our two-dimensional map. This general form of the matrix is 



Df 



df a /dm a df u /dm T 
df T /dm a df T /dm T 

Pa J (sech 2 % (jm a + r)^ \ E s =±i t&nh ft a ( J m a s + I] 
\ E s =±i tanh h{Jm T s + 1) P T j(sech 2 (3 T (jm T + a J ^ 

For the GP phase, at m a = m T = 0, the Jacobian matrix simplifies to 



(25) 



Df 



with eigenvalues 



where 



P a Jsech 2 /3 CT tanh/3 CT 
tanh/3 r /3 T Jsech 2 /9 r 



X ± = 1 -(-b±VA), 



J [ P a sech 2 P a + /3 r sech 2 /3. 



T I ) 

2 



A = J p a sech z p a - p r sech z p T + 4 tanh p a tanh fi T . 



(26) 

(27) 

(28) 
(29) 



For P T = Pa = 0, both eigenvalues are zero indicating that the GP phase is a sink, i.e., 
stable to perturbations. For small, but non-zero inverse temperatures (high-temperature 
limit), we have 



X ± -i:[J(Pa + Pr)±\J(Pa-Pr 



(30) 



which gives eigenvalues J P a and jp r . As stability in the direction of each of the eigenvec- 
tors depends on the magnitude of the corresponding eigenvalue, we see that the GP phase 
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FIG. 5. Stability of the GP phase for high temperatures (small /?'s). 



becomes unstable below some temperature value. For either (3 a > l/\J\ or f3 T > 1/\J\ sepa- 
rately, it becomes a saddle. If both conditions are satisfied it becomes a source. Note that 
these conditions are exactly the ones defining the critical point of the naive mean field Ising 
model. The results for high temperatures are indicated in Fig [5j 
The results for low temperatures with (3 a , (3 T — > oo are 



and the GP is always a saddle. 

An interesting case occurs in the regime where one of the temperatures is high and the 
other is low. By taking /3 r — > and /3 a — > oo, straightforward calculations show that 



For sufficiently large values this points to a stability of the GP phase for all interaction 
values. This can be seen from the basin of attraction depicted in Fig. [6j The figure shows 
the basin of attraction for fixed points of the two-dimensional map f[T3"j) with parameters 
given by J = 2, J' = 1, /3 a = 10 3 and /3 T = 0. The three black circles correspond to the 
stable fixed points and the shaded regions around each of them are their respective basins 
of attraction. There are also two unstable fixed points, located in the same horizontal line 
as the stable ones at the border of the touching basins of attraction. 



X ± « ±1 + 2J(/3 CT e- 2 ^ + /3 r e" 2 ^), 



(31) 



X ± w J{2f3 a e~ 2 ^ + /3 T ) ± t/WIFt. 



(32) 
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FIG. 6. Basin of attraction for values J = 2, J' = 1, f} a = 10 3 and f3 T = 0. 



V. CONCLUSIONS 



With the aim of understanding the behaviour of non-equilibrium systems with two- 
temperatures, we introduced and solved a simple model of a system comprising two spin- 1/2 
components, referred to as the a and r components; these are coupled to separate thermal 
reservoirs held at different temperatures T a and T T , respectively. For visualisation purposes, 
we considered the two system components as belonging to two different planes. Each com- 
ponent has then an intraplane long-ranged mean-field interaction J while it interacts with 
the other component via an interplane local interaction J'. 

Using the generating functional technique we have found a pair of exact, self-consistent 
coupled equations for the magnetisations of the two components m" and m r , the solution 
of which is given by the fixed-points of the resulting two-dimensional non-linear map. This 
map depends on three degrees of freedom which are given by the ratio between the inter and 
intra-component couplings J — J/ J' and the scaled inverse temperatures of each component 
a = J'(3 a and T = J'fi T . 

The equations for this model are reminiscent of equilibrium equations for mean-field 
magnetic systems and allow for numerical solution of the fixed-points at any parameter 
values. Analytically, we have been able to study the possible phases of this system in the 
limits of low and high temperature, where one could identify paramagnetic and ferromagnetic 
phases, showing that these are well defined steady states. 
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The analysis of the particular phase where both components have zero magnetisation, 
which we termed the Global Paramagnetic (GP) phase, show that the stability of this phase 
has some interesting characteristics, especially in the low and high temperature limits. We 
have not found any evidence for chaotic behaviour in general, but have no definitive proof 
either for its existence or absence as each orbit would need to be analysed individually. We 
conjecture that this is due to the nature of the stochastic dynamics that give rise to this 
mapping as the magnetisations fluctuate very close to the mean values at the large system 
limit which we investigate. 

The framework introduced and the results obtained can be easily modified to accommo- 
date other cases; one obvious modification being the spatial nature (local/non-local) of the 
interactions. Also, with the appropriate changes, one can apply this theory to several other 
non-equilibrium practical situations where a steady state is reached. 

One of these cases, which we are currently exploring, is related to astrophysical dark 
matter as the dark sector might have, in additional to its gravitational coupling to the 
visible sector, exclusive dark matter interactions. This opens up the possibility of systems at 
different temperatures to co-exist at the same physical location but in different sectors. These 
systems can then be heated independently by processes on their own sectors and exchange 
energy only via gravitational interactions. The current framework could be modified to 
accommodate long-range gravitational "interplane" interaction and short-rang "intraplane" 



interaction according to bounds obtained from astrophysical observations |2lj. Analysis of 
the phases and the heat-transfer between both sectors could shed light on the nature of these 
potential interactions or serve as an additional detecting mechanism. 

Another extension of this framework, which is also under way, is the application to 



superlattices 23[. Superlattices are metamaterials constructed of alternate layers of two 
different materials that hold useful properties. In a forthcoming study, we appropriately 
alternate thermal conducting magnetic layers with thermal isolating ones such that each 
one of n > 2 conducting layers can be coupled to a different thermal reservoir. This would 
result in an n-temperatures systems and the magnetisations would then be described by an 
n-dimensional non-linear map with the promise of a rich and interesting phase diagram. 

We believe that there are potentially many cases of non-equilibrium physical systems 
comprising a number of components, each of which is held at a different temperature and 
exhibits interactions of a different nature to the cross-component ones. This analysis of the 
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particularly simple case of interacting magnetic systems held at different temperatures, can 
be carried out rigorously and paves the way to better understanding of similar cases in a 
variety of fields. 
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